function [ output, roiX, roiZ ] = tdSAFT2D( scanData, varargin )
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
%
scanDataValidationFcn = @(x) isstruct(x) && isfield(x, 'CscanData') && ...
    isfield(x, 'timeScale') && isfield(x, 'AscanValues');
validOutput = {'dB','normalized','raw'};
checkOutput = @(x) any(validatestring(x,validOutput));
validPostproc = {'env','rectified','none'};
checkPostproc = @(x) any(validatestring(x,validPostproc));
validApodization = {'none','rectangular','Hann'};
checkApodization = @(x) any(validatestring(x,validApodization));

% Verify scanData parameter
assert(scanDataValidationFcn(scanData), 'scanData invalid structure');

% Get scan data informations
Tresol = (scanData.timeScale(2) - scanData.timeScale(1))*1e-6;
Tini = scanData.timeScale(1)*1e-6;
Tmax = scanData.timeScale(end)*1e-6;
if strcmp(scanData.CscanData.WaveType, 'L')
    c = scanData.CscanData.Cl;
else
    c = scanData.CscanData.Cs;
end
ERMv = c/2;
iERMv = 2/c;
maxVal = scanData.CscanData.MaxSignal;
scanX = scanData.CscanData.X*1e-3;
scanY = scanData.CscanData.Y*1e-3;
tamX = length(scanX);
tamY = length(scanY);
deltaX = scanX(2) - scanX(1);
minZ = 0;
maxZ = scanData.timeScale(end)*1e-6*ERMv;
deltaZ = Tresol*ERMv;
apertAngle = asin(0.51*c/(scanData.CscanData.Frequency*1e6*scanData.CscanData.ProbeDiameter*1e-3));
filter=true;
% Create parser to process arguments
inpP = inputParser;
% ROI X parameters
addParamValue(inpP, 'minX', scanX(1), @isnumeric);
addParamValue(inpP, 'maxX', scanX(end), @isnumeric);
% ROI Z parameters
addParamValue(inpP, 'minZ', minZ, @isnumeric);
addParamValue(inpP, 'maxZ', maxZ, @isnumeric);
% Plot output
addParamValue(inpP, 'plotOut', true, @islogical);
% Debug output
addParamValue(inpP, 'debug', false, @islogical);
% Y-slice selector
addParamValue(inpP, 'Yslice', find(scanY==0), @isnumeric);
% Output style
addParamValue(inpP, 'output', validOutput{1}, checkOutput); 
% dB limit
addParamValue(inpP, 'dBlim', -26, @isnumeric);
% X step
addParamValue(inpP, 'stepX', 1, @isnumeric);
% Post process
addParamValue(inpP, 'Postproc', validPostproc{1}, checkPostproc);
% Profile plot
addParamValue(inpP, 'profile', true, @islogical);
% Apodization function
addParamValue(inpP, 'apodization', validApodization{1}, checkApodization);
% Apodization window length
addParamValue(inpP, 'apod_length', 1.0, @isnumeric);

% Parse input arguments
parse(inpP, varargin{:});

% Define roiX
roiX = inpP.Results.minX:deltaX:inpP.Results.maxX;
tamRoiX = length(roiX);

% Define roiZ
if(isfield(scanData.CscanData, 'DepthBW'))
    roiZ = inpP.Results.minZ:deltaZ:scanData.CscanData.DepthBW*1e-3;
else
    roiZ = inpP.Results.minZ:deltaZ:inpP.Results.maxZ;
end
tamRoiZ = length(roiZ);

% Verify Y-slice
if(tamY > 1)
    assert(inpP.Results.Yslice > 0 && inpP.Results.Yslice <= tamY, ...
        'Yslice parameter out of range');
    yi = inpP.Results.Yslice;
else
    yi = 1;
end

% Get X step
stepX = inpP.Results.stepX;
stepMsg = sprintf('X step: %d', stepX);

output = zeros(tamRoiZ, tamRoiX);
if(inpP.Results.debug)
    tic;
end

% Set Stolt migration parameters
[nt0,nx0] = size(scanData.AscanValues);

% Zero-padding
nt = 2^(nextpow2(nt0)+1);
nx = 2*nx0;

% Get dt and dx
dt = (scanData.timeScale(2) - scanData.timeScale(1))*1e-6;
dx = (scanData.CscanData.X(2) - scanData.CscanData.X(1))*1e-3;

if(inpP.Results.debug)
    tic;
end

% Calculate the matched filter
d = scanData.CscanData.ProbeDiameter*1e-3;
fc = scanData.CscanData.Frequency*1e6;
bw = scanData.CscanData.Bandwidth*1e6;
[A, Hed] = calcMatchedFilter(nt, 1/dt, nx, dx, d, fc, bw);

% Calculate FFT of scanned data
fftRF = fftshift(fft2(scanData.AscanValues, nt, nx));

% Apply matched filter
if(inpP.Results.filters)
    fftRF = bsxfun(@times, fftRF, Hed');
    ep = 1e1;
    AA = A./(A.^2 + ep^2);
    fftRF = bsxfun(@times, fftRF, AA);
end

% Calculate Stolt migration
imgRF = stoltmig(fftRF, 1/dt, dx, ERMv, scanData.CscanData.BeamAngle*pi/180);

% IFFT & Migrated RF
migRF = ifft2(ifftshift(imgRF),'symmetric');
output = migRF(1:nt0,1:nx0);

if(inpP.Results.debug)
    toc
end

% Output post-process
switch inpP.Results.Postproc
    case 'env'
        output = abs(hilbert(output));
        envMsg = 'w-k SAFT-2D envelope image';
        profMsg = 'w-k SAFT-2D envelope profile';
    case 'none'
        envMsg = 'w-k SAFT-2D image';
        profMsg = 'w-k SAFT-2D profile';
end

switch inpP.Results.output
    case 'raw'
        titleProfMsg = {[profMsg ' (raw)'], ...
            sprintf('Y Slice: %d', inpP.Results.Yslice)};
        titleEnvMsg = {[envMsg ' (raw)'], ...
            sprintf('Y Slice: %d', inpP.Results.Yslice)};
        imgLimits = [min(min(output)) max(max(output))];
    case 'dB'
        output = 20*log10(abs(output) ./ max(max(abs(output))));
        titleProfMsg = {[profMsg sprintf(' (%.2f dB)', inpP.Results.dBlim)], ...
            sprintf('Y Slice: %d', inpP.Results.Yslice)};
        titleEnvMsg = {[envMsg sprintf(' (%.2f dB)', inpP.Results.dBlim)], ...
            sprintf('Y Slice: %d', inpP.Results.Yslice)};
        imgLimits = [inpP.Results.dBlim 0];
    case 'normalized'
        output = output ./ max(max(abs(output)));
        titleProfMsg = {[profMsg ' (normalized)'], ...
            sprintf('Y Slice: %d', inpP.Results.Yslice)};
        titleEnvMsg = {[envMsg ' (normalized)'], ...
            sprintf('Y Slice: %d', inpP.Results.Yslice)};
        imgLimits = [0 1];
end

% Plot results
if(inpP.Results.plotOut)
    if(inpP.Results.profile)
        figure; plot(roiX, max(output));
        title(titleProfMsg);
        xlabel('X (m)');
        ylabel('max value');
        grid on;
    end
    figure; imagesc(roiX, roiZ, output, imgLimits);
    colormap(jet);
    colorbar;
    title(titleEnvMsg);
    xlabel('X (m)');
    ylabel('Z - depth (m)');
    grid on;
end
end
